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The image of planetary nebulae (PN), supernova remnant (SNR) and Eta-Carinae is made 
by three different physical processes. The first process is the expansion of the shell that can 
be modeled by the canonical laws of motion in the spherical case and by the momentum 
conservation when gradients of density are present in the interstellar medium. The quality 
of the simulations is introduced along one direction as well along many directions. The 
second process is the diffusion of particles that radiate from the advancing layer. The 3D 
diffusion from a sphere , the ID diffusion with drift and ID random walk are analyzed. The 
third process is the composition of the image through an integral operation along the line 
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of sight. The developed framework is applied to three PN which are A39 , the Ring nebula 
and the etched hourglass nebula MyCn 18, the hybrid object Eta-Carinae , and to two SNR 
which are SN 1993J and SN 1006. In all the considered cases a careful comparison between 
the observed and theoretical profiles in intensity is done. 



1. Introduction 

The spherical explosion in galactic astrophysics models two different objects: 

• The planetary nebula (PN) that are characterized by small terminal velocities of the 
order of few km/s and small energy , E, released in the expansion , E w 10^^ erg. 

• The supernova remnant (SNR) that are characterized by high velocities 5000A;m/s 
and high energy involved , E k, IQ^^ erg. 

These two main classifications does not cover peculiar astrophysical objects such as the 
nebula around 77-Carinae which is characterized by high velocity , f« 300 km/s and low 
energy involved , E 10'^^ erg. We now summarize the existing models on these three 
astrophysical objects. PN The PN rarely presents a circular shape generally thought to be 
the projection of a sphere on the sky. In order to explain the properties of PN, [1] proposed 
the interacting stellar wind (ISW) theory. Later on \2j proposed the two wind model and the 
two phase model. More often various types of shapes such as elliptical , bipolar or cigar are 
present, see |[3]-[8|. The bipolar PNs , for example , are explained by the interaction of the 
winds which originate from the central star , see [9-12|. Another class of models explains 
some basic structures in PNs through hydrodynamical models, see |,13„14] or through self- 



organized magnetohydrodynamic (MHD) plasma configurations with radial flow, see 1 15 1. 
An attempt to make a catalog of Une profiles using various shapes observed in real PNs was 
done by |[16|. This ONLINE atlas , available at 

http://132.248. 1.102/Atlas_profiles/img/, is composed of 26 photo-ionization models corre- 
sponding to 5 geometries, 3 angular density laws and 2 cavity sizes, four velocity fields for 



a total of 104 PNs, each of which can be observed from 3 different directions. 1 17 1 suggest 
that a planetary nebula is formed and evolves by the interaction of a fast wind from a central 
star with a slow wind from its progenitor , an Asymptotic Giant Branch (AGB) star. 



ry-Carinae The nebula around r/-Carinae was discovered by |18| and the name "the 



Homunculus" arises from the fact that on the photographic plates it resembled a small 



plump man, see 1 19 1. More details on the various aspects of r^-Carinae can be found in [20]. 
The structure of the Homunculus Nebula around ?7-Carinae has been analyzed with different 
models, we cite some of them: 

• The shape and kinematics is explained by the interaction of the winds expelled by the 
central star at different injection velocities, see [9,]. 

• The possibility that the nebulae around luminous blue variables (LBVs) are shaped 



by interacting winds has been analyzed by |21 1. In this case a density contrast profile 
of the form p = pq{\ + b cos^ 0) was used where is the angle to the equatorial 
plane. 
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The origin and evolution of the bipolar nebula has been modeled by a numerical 
two-dimensional gasdynamic model where a stellar wind interacts with an aspherical 
circumstellar environment, see [10]. 

Cooling models form ballistic flows (that is, a pair of cones each with a spherical 
base) whose lateral edges become wrinkled by shear instabilities, see |[22|. 



The scaling relations derived from the theory of radiatively driven winds can model 
the outflows from luminous blue variable (LBV) stars, taking account of stellar rota- 
tion and the associated latitudinal variation of the stellar flux due to gravity darken- 
ing. In particular for a star rotating close to its critical speed, the decrease in effective 
gravity near the equator and the associated decrease in the equatorial wind speed re- 
sults naturally in a bipolar, prolate interaction front, and therefore in an asymmetric 
wind, see p3[ . 

Two oppositely ejected jets inflate two lobes (or bubbles) representing a unified model 
for the formation of bipolar lobes, see |24j25). 



A two-dimensional, time-dependent hydrodynamical simulation of radiative cooling, 
see 12611. 



Launch of material normal to the surface of the oblate rotating star with an initial 
kick velocity that scales approximately with the local escape speed, see ||27|. 



A 3D model of wind-wind collision for X-ray emission from a supermassive star, 
see |j28J. 

Two-dimensional hydrodynamical simulations of the eruptive events of the 1840s 
(the great outburst) and 1890s (the minor outburst), see p9|. 



SNR The study of the supernova remnant (SNR) started with 1 30 1 where an on ongo- 
ing coUisional excitation as a result of a post-explosion expansion of the SNR against the 
ambient medium was suggested. The next six decades where dedicated to the deduction 
of an analytical or numerical law of expansion. The target is a relationship for the instan- 
taneous radius of expansion, R, of the type oc where t is time and m is a parameter 
that depends on the chosen model. On adopting this point of view, the Sedov expansion 
predicts R oc t^'"^, see [31], and the thin layer approximation in the presence of a constant 



density medium predicts R cx t^-"^^, see [|32[. A simple approach to the SNR evolution in 



the first 10"^ yr assumes an initial free expansion in which R oc t until the surrounding mass 
is of the order of 1 Mq and a second phase characterized by the energy conservation in 
which according to the Sedov solution R oc t^/^, see {33} . A third phase characterized by 



an adiabatic expansion with R oc t^/^ starts after 10^ yr, see [33]. A more sophisticated 



approach given by [ [34| 35 1 analyzes self-similar solutions with varying inverse power law 



exponents for the density profile of the advancing matter, R and ambient medium, R 



n-3 



The previous assumptions give a law of motion R oc t"-^ when n > 5. Another example 



is an analytical solution suggested by ] 36 ] where the radius-time relationship is regulated 
by the decrease in density: as an example, a density proportional to R~^ gives R oc t^/^. 
With regard to observations, the radius-time relationship was clarified when a decade of 
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very-long-baseline interferometry (VLBI) observations of SN 1993J at wavelengths of 3.6, 
6, and 18 cm became available, see |j37|-[39|. As a first example, these observations col- 
lected over a 10 year period can be approximated by a power law dependence of the type 
R oc t^'^"^. This observational fact rules out the Sedov model and the momentum con- 
servation model. In this paper we describe in Section [T] the observed morphologies of 
PNs, SNRs and 7] — car. Section |3]] analyzes five different laws of motion that model the 
spherical and aspherical expansion and Section [4j] applies the various law of motion to well 
defined astrophysical objects introducing the quality of the simulation. Section [5T] reviews 
old and new formulas on diffusion , Section [6] reviews the existing situation with the ra- 
diative transport equation and Section [T] contains detailed information on how to build an 
image of the astrophysical objects here considered. 



2. Astrophysical Objects 

This section presents the astronomical data of a nearly spherical PN known as A39, a 
weakly asymmetric PN , the Ring nebula , and a bipolar PN which is the etched hourglass 
nebula MyCn 18 . The basic data of ?7-Carinae which is not classified as a PN due to the 
high velocities observed are also reported. The section ends with two SNR , the symmetric 
SN 1993 J and the weakly asymmetric SN 1006 . 



2.1. A circular PN , A39 

The PN A39 is extremely round and therefore can be considered an example of spherical 



symmetry, see for example Figure 1 in |40 1 . In A39 the radius of the shell , Rsheii is 

Rsheii = 2 A2 X 10^^ Qrr D21 cm = 0.78 PC , (1) 



where 677 is the angular radius in units of 77" and D21 the distance in units of 2.1 kpc , 
see |40| . The expansion velocity has a range [32 o 37 — ] according to |41 1 and the age 



of the free expansion is 23000 yr, see |40 1. The angular thickness of the shell is 



6rsheii = 3.17 10i^eio£'2icm = 0.103 pc 



(2) 



where ©10 is the thickness in units of 10.1" and the height above the galactic plane is 
1.42 kpc , see [40]. The radial distribution of the intensity in [OIII] image of A39 after 
subtracting the contribution of the central star is well described by a spherical shell with a 
10" rim thickness, see Figure [T]and ||40j. 



2.2. A weakly asymmetric PN , M57 

The Ring nebula , also known as M57 or NGC6720 , presents an elliptical shape charac- 
terized by a semi-major axis of 42", a semi-minor axis of 29.4" and ellipticity of 0.7, see 



Table I in [42]. The distance of the Ring nebula is not very well known ; according to |43 1 
the distance is 705 pc . In physical units the two radii are 

Rshell,minoT = 0. 1029.4-D7O5 pc scTJii - minor radius 
Rsheii,major = O.14042-D7O5 pc scmi - major radius , (3) 
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Arcsec from Nebula Center 



Figure 1. Cut of the relative intensity of the PN A39 crossing the center in the east-west 
direction (dotted line with some error bar) and the rim model (full line) fully described in 
Jacoby et al. (2001) . 

where 029.4 is the angular minor radius in units of 29.4", 642 is the angular major radius 
in units of 42" and -D705 the distance in units of 705 pc. The radial velocity structure in 
the Ring Nebula was derived from observations of the H2 (molecular Hydrogen) v = 1- 
S(l) emission line at 2.122 fim obtained by using a cooled Fabry- Perot etalon and a 
near-infrared imaging detector , see |42| . The velocity structure of the Ring Nebula covers 
the range [-30.3 o 48.8 ^] . 



2.3. A strongly asymmetric PN , MyCn 18 

MyCn 18 is a PN at a distance of 2.4 kpc and clearly shows an hourglass-shaped nebula, 
see [ i44|45J . On referring to Table 1 in | |46| we can fix the equatorial radius in 2.80 x 10^^ cm 
, or 0.09 pc , and the radius at 60° from the equatorial plane 3.16 x 10^^ cm or 0.102 pc . 
The determination of the observed field of velocity of MyCn 1 8 varies from an overall value 



of 10 ^ as suggested by the expansion of [01 1 1] , see |45| , to a theoretical model by |46 
in which the velocity is 9.6 ^ when the latitude is ° (equatorial plane) to 40.9 ^ when 
the latitude is 60 °. 



2.4. Homunculus nebula 

The star 77-Carinae had a great outburst in 1840 and at the moment of writing presents a 
bipolar shape called the Homunculus, its distance is 2250 pc, see [47 J . A more refined clas- 
sification distinguishes between the large and little Homunculus, see HSl. The Homunculus 



has been observed at different wavelengths such as the ultraviolet and infrared by |20 49 1, 
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Time (yr) 



Figure 2. Radius in pc versus year of the SNR SN 1993J with vertical error bars. The em 
bands are A = 3.6 cm and A = 6 cm. The data are extracted from Table 1 in Marcaide et 
al. 2009. The dotted line represents an expansion at a constant velocity. 



x-ray by fSO^, [Fe//]A16435 by fST), ammonia by (52^, radio-continuum by |48|, near- 
infrared by [53] and scandium and chromium lines by |54|. 

Referring to Table 1 in [55 1, we can fix the major radius at 22014 AU (0.106 pc) and 
the equatorial radius at 2100 AU (0.01 pc). The expansion velocity rises from 93 km/ s 
at the equator to 648 km /s in the polar direction, see Table 1 and Figure 4 in f55l. The 
thickness of the H2 shell is roughly 2 — 3% of the polar radius, see ||55J. 

2.5. A circular SNR, SN 1993J 



The supernova SN1993J started to be visible in M81 in 1993, see |56|, and presented a 



see |57|. The expansion of SN 1993 J has been monitored in various bands over a decade 



circular symmetry for 4000 days, see |37|. Its distance is 3.63 Mpc (the same as M81), 
see |57|. The expansion of SN 1993J has 1 
and Figure |2] reports its temporal evolution. 
Dser 

V = 8474 ^ at t = 10.53 yr. We briefly recall that |58 1 quote an inner velocity from the 



The observed instantaneous velocity decreases from v = 15437 ^ at t = 0.052 yr to 



shapes of the lines of « 7000^ and an outer velocity of k. 10000 
2.6. A weakly asymmetric SNR, SN 1006 

The diameter of the known remnants spans the range from 3 pc to 60 pc and attention is fixed 
on SN 1006 which started to be visible in 1006 AD and its possible diameter of 12.7 pc. 



km 
s 



see p9| . More precisely, on referring to the radio-map of SN 1006 at 1370 MHz by | |60] , 
it can be observed that the radius is greatest in the north-east direction. From the radio- 
map previously mentioned we can extract the following observed radii, R = 6.8 pc in the 
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polar direction and R = 5.89 pc in the equatorial direction. Information on the thickness of 
emitting layers is contained in a recent study by |61 1 which analyzes Chandra observations 
(i.e., synchrotron X-rays) from SN 1006 . The observations found that sources of non- 
thermal radiation are likely to be thin sheets with a thickness of about 0.04 pc upstream and 
0.2 pc downstream of the surface of maximum emission, which coincide with the locations 
of Balmer-line optical emission , see |62| . The high resolution XMM -Newton Reflection 
Grating Spectrometer (RGS) spectrum of SN 1006 gives two solutions for the O VII triplet. 
One gives a shell velocity of 6500 km/s and the second one a shell velocity of 9500 
km/ s, see |j63j|. The value here adopted for the magnetic field can he. H = 10^ Gauss as 
suggested by ||6T|. 



3. Law of motion 



This Section presents three solutions for the law of motion that describe a symmetric ex- 
pansion. The momentum conservation is then applied in cases where the density of the 
interstellar medium is not constant but regulated by exponential or power law behavior. 



3.1. Spherical Symmetry - Power law solution 

The equation of the expansion of the SNR can be modeled by a power law of the type 

R{t)=Ro{^r , (4) 

where R is the radius of the expansion, t is the time , Rq is the radius at i = to and a an 
exponent which can be found from the numerical analysis. In the case of SN 1993 J we have 
a =0.828 , i?o=0.0087 pc and to = 0.498 yr 



3.2. Spherical Symmetry - Sedov solution 

The momentum conservation is applied to a conical section of radius R with a solid angle 
A ri, in polar coordinates, see | [33| 

j^{AMR) = AF, (5) 

where 

rR 

AM= / p{R,6,4))dV , (6) 
Jo 

is the mass of swept-up interstellar medium in the solid angle A il, p the density of the 
medium , P the interior pressure and the driving force: 

AF = PR^ An . (7) 

After some algebra the Sedov solution is obtained, see 1 3 1 33] 

2b Et'-^ ^" 



B(«)=|^ — I . (8) 
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where E is the energy injected in the process and t the time. 
Another slightly different solution is formula (7.56) in 132 



where the difference is due to the adopted approximations. 

Our astrophysical units are: time (t^), which is expressed in 10^ yr units; £'42, the 
energy in W^^ erg; and no the number density expressed in particles cm^^ (density p = 
nom, where m=1.4mH). With these units equation ([8]) becomes 

/ p ^ 2\ 1/5 

0.198 pc . (10) 



The expansion velocity is 
which expressed in astrophysical units is 



no 



V(*)-"46-^^ . (12) 



Equations ( 10 1 and ( 12 1 represent a system of two equations in two unknowns : and i?42 



. By inserting for example i? = 0.78 pc in equation ( 10) we find 



t4 = 77.15^^ , (13) 



and inserting V = 35 km s ^ in equation ( 12 1 we obtain 



0.3954 Y^E^ = 35 . (14) 



The previous equation is solved for £'42 = 7833.4 that according to equation ( [13] ) means 
t4=.87173. These two parameters allows a rough evaluation of the mechanical luminosity 
L = y that turns out to be L 2.847 10^^^ergs s~^. This value should be bigger than the 
observed luminosities in the various bands. As an example the X-ray luminosity of PNs , 
Lx, in the wavelength band 5-28 A has a range [10^*^-^ ^ s~^] , see Table 3 

in||64). 

Due to the fact that is difficult to compute the volume in an asymmetric expansion the 
Sedov solution is adopted only in this paragraph. 

3.3. Spherical Symmetry - Momentum Conservation 

The thin layer approximation assumes that all the swept-up gas accumulates infinitely in a 
thin shell just after the shock front. The conservation of the radial momentum requires that 

^7tR^pR = Mo , (15) 



Interaction of PN and SNR with the ISM 



9 



where R and R are the radius and the velocity of the advancing shock , p the density of the 
ambient medium , Mq the momentum evaluated att = , Rq the initial radius and Rq the 



initial velocity , see |32 65 1. The law of motion is 



R = Roll + A^^{t-to) 



(16) 



and the velocity 



R = Ro (l + 4|^(*-*o) 



(17) 



From equation ( 16 1 we can extract Rq and insert it in equation ( 17 1 



R 



It — JXn 



4(t - to) Rl 



1 + 



It — Itn 



Rq 



(18) 



The astrophysical units are: and to,4 which are t and expressed in 10^ yr units, Rpc and 
i?o,pc which are R and Rq expressed in pc, Rkms and Ro,kms which are R and i?o expressed 
in ^ . Therefore the previous formula becomes 



_ 3 

1 p4 _ r4 / r4 _ p4 \ 4 



(19) 



0,pc 



0,pc 



On introducing i?o,pc = 0.1 , Rpc = 0.78 , Rkms = 34.5^ , the approximated age of A39 
is found to be t4 — to,4 = 50 and Ro^kms = 181.2. 



3.4. Asymmetry - Exponential medium 

Given the Cartesian coordinate system (x, y, z) , the plane z = will be called equatorial 
plane and in polar coordinates z = R sin(^), where 9 is the polar angle and R the distance 
from the origin . The presence of a non homogeneous medium in which the expansion takes 
place can be modeled assuming an exponential behavior for the number of particles of the 
type 

z R X sm(9) 
n(z)=noexp— — =noexp , (20) 

where R is the radius of the shell , no is the number of particles at i? = i?o and h the scale. 
The 3D expansion will be characterized by the following properties 

• Dependence of the momentary radius of the shell on the polar angle 9 that has a range 

[-90° o +90°]. 

• Independence of the momentary radius of the shell from cj) , the azimuthal angle in 
the x-y plane, that has a range [0° o 360°] . 
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The mass swept, M, along the solid angle A 17, between and R is 

M{R) = ^munoImiR) + ^^Rlnomu , (21) 

where 

I^{R)= [''r'e^v-'-^dr , (22) 
where Rq is the initial radius and mn the mass of the hydrogen . The integral is 

h{2h'^ + 2 Rohsin (6) + Ro"^ (sin (0))^) ^.teW 

h(2h'^ + 2Rh sin (9) + R^ (sin {9)f) e~ 

^ 5 • (23) 

(sin(0))^ 

The conservation of the momentum gives 

M{R)R = M{Ro)Ro , (24) 

where R is the velocity at R and Rq the initial velocity at R = Rq. 

In this differential equation of the first order in R the variable can be separated and the 
integration term by term gives 

/ M{r)dr = M{Ro)Ro x {t - to) , (25) 

where t is the time and to the time at Rq. The resulting non linear equation Tnl expressed 
in astrophysical units is 

_ Ro,pc sin(9) _ Ro,pc Bin(9) 

J'NL = -Qe V hpc^-hpce V {sin {9) f Ro,pc^ 

_ R0,pc sin(9) _ Rp^pc Bin(9) 

-6V^e V sin (9) Ro,pc- 3 hpc'^e V {sin {9)y Rg^pc"^ 

_ Rpc sin(fl) _ Rpc Bin(9) 

- Ro, pepsin {6)f + Qe V v4 + 4e V /i^.^^^^ (0) 

+e V /ip,2^j,,2(sin(^))2 

_ Ro,pc sin(f)) ^ flp.pc Bin(9) 

+2e V sin(e) +2e V hp^^ Rp^ {sin {e)Y Ro,pc 



_ Ro,pc sin(9) 

+e V hpc Rpc {sin {9) f Ro,pc'^ 



+ (sin (0))^ Ro,pc^Rpc - 0.01 (sin (^))^ Ro,pc^Ro,kms {U - to,4) = , (26) 

where ^4 and to,4 are t and to expressed in 10^ yr units, Rpc and i?o,pc are and Rq 
expressed in pc, Rkms and Ro,knis are ii and Rq expressed in ^,9 is expressed in radians 
and hpc is the the scale , h , expressed in pc. It is not possible to find Rpc analytically and a 
numerical method should be implemented. In our case in order to find the root of Tnl, the 
FORTRAN SUBROUTINE ZRIDDR from |66| has been used. 

The unknown parameter — to,4 can be found from different runs of the code once 
Ro^pc is fixed as w 1/10 of the observed equatorial radius , Ro,kms is 200 or less and hpc ~ 
2 X Ro,pc- 
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3.5. Asymmetry - Power law medium 

A possible form for a power law profile of the medium surrounding the Homunculus nebula 
is 

71(2;) = no(— , (27) 

where z = R x sin(^?) is the distance from the equatorial plane, R is the instantaneous 
radius of expansion, no is the number of particles at R = Rq, Ro is the scale and q is a 
coefficient > 0. 

The swept mass, M, along the solid angle A Q between and R is 

M{R) = ^lAmHnoIm{R) + ^TTRlnomHlA , (28) 

where 

ImiR)= Tr^-^rdr , (29) 
JRo ^0 

where Rq is the initial radius and ttih is the mass of hydrogen. The integral is 

Im{R)= \ ^ . (30) 

6 — a 

Conservation of momentum gives 

M{R)R = M{Ro)Ro , (31) 

where R is the velocity at R and Rq is the initial velocity at R = Rq , M{R) and M(Rq) 
are the swept masses at R and Rq respectively 

In this first-order differential equation in R, the variables can be separated. Integration 
term-by-term gives 

/ M{r)dr = M{Ro)Ro x {t - to) , (32) 

where t is the time and to is the time at Rq. The resulting non-linear equation Tnl ex- 
pressed in astrophysical units is 

Jtvl = - Ro,pc^<y^ + Rpc Ro,pc^Oi^ + Rpc Ro,pc^ (sin {9))~ " a 
+7.0 Ro,pc^a - Ro,pc^ (sin {6))- " a - 7.0 R^c Ro,pc^a 
+3 i?o,pc' (sin (0))- " - 12 Ro,pc^ - 4.0 Rp^ Ro,pc^ (sin {9)^ 
+ Rpc ""'"^■''-Ro.pc" (sin (0)) + 12 Rpc Ro,pc'^ 
— 0.122 Ro^pc^Ro,kms {ti — to^4) — 0.01 Ro,pc^Ro,kms {ti — to,4)ot^ 

+0.0714 Ro,pc^Ro,kms {U - to,4 )a = , (33) 

where and to,4 are t and to expressed in 10^ yr units, Rpc and i?o,pc are R and Rq 
expressed in pc, Rkms and RQ^kms are R and Rq expressed in ^ and 9 is expressed in 
radians. 

It is not possible to find Rpc analytically and a numerical method must be implemented. 
In our case, in order to find the root of Tnl, the FORTRAN SUBROUTINE ZRIDDR 
from [66] has been used. The unknown parameters, -Ro,pc and Ro^kms^ are found from 
different runs of the code, — to,4 is an input parameter. 
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Table 1 . Data of the simulation of the PN Ring nebula 



Initial expansion velocity ,Ro,kms 200 

Age(t4-to,4) ' 0.12 

Initial radius Ro,pc 0.035 

scaling hpc 2 x Rq , 



4. Applications of the law of motion 

From a practical point of view, e , the percentage of reliability of our code can also be 
introduced, 

e = (1 - l(-^PC,obs-i?pc,num)| ^ _ ^ ^34^ 
Rpc,ohs 

where Rpc,ohs is the radius as given by the astronomical observations in parsec , and Rpc,num 
the radius obtained from our simulation in parsec. 

In order to test the simulation over different angles, an observational percentage of 
reliability ,eobs> is introduced which uses both the size and the shape, 

.r,nr-i Ylj\^pc,ohs - Rpc,num\j 

eobs = 100(1 ^ ), (35) 

/ y j -^pCjOhs j 

where the index j varies from 1 to the number of available observations. 
4.1. Simulation of a PN , Ring nebula 

A typical set of parameters that allows us to simulate the Ring nebula is reported in Table[T] 
The complex 3D behavior of the advancing Ring nebula is reported in Figure [3] and 
Figure |4] reports the asymmetric expansion in a section crossing the center. In order to 
better visualize the asymmetries Figure [5] and Figure [6] report the radius and the velocity as 
a function of the position angle 6. The combined effect of spatial asymmetry and field of 
velocity are reported in Figure |7] 



The efficiency of our code in reproducing the observed radii as given by formula ( 34 ) 
and the efficiency when the age is five time greater are reported in Table |2] 

An analogous formula allows us to compute the efficiency in the computation of the 
maximum velocity , see Table [3] 

4.2. Simulation of PN , MyCn 18 

A typical set of parameters that allows us to simulate MyCn 18 is reported in Table [4] 

The bipolar behavior of the advancing MyCn 18 is reported in Figure [8] and Figure |9] 
reports the expansion in a section crossing the center. It is interesting to point out the simi- 
larities between our Figure |9]of MyCn 18 and Figure 1 in [16J which define the parameters 
a and h of the Atlas of synthetic line profiles. In order to better visualize the two lobes 



Figure 10 reports the radius as a function of the position angle 9. 
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Figure 3. Continuous three-dimensional surface of the PN Ring nebula : the three Eule- 
rian angles characterizing the point of view are $=180 °, 0=90 ° and ^'=-30 °. Physical 
parameters as in Table [T] 



Table 2. Reliability of the radii of the PN Ring nebula. 





i?up(pc) polar direction 


i?eq(pc) equatorial plane 


^obs 


0.14 


0.1 


-Rnum(our code) 


0.125 


0.102 




89 


97 


e (%) for a time 5 times greater 


27 


41 



14 



Zaninetti 



IN 



' 



-ai -005 



X [pc] 



006 ai 



Figure 4. Section of the PN Ring nebula on the x-z. plane. The horizontal and vertical axis 
are in pc. Physical parameters as in Table [T] 



Table 3. Reliability of the velocity of the PN Ring nebula 



y( — ) maximum velocity 



Kbs 48.79 

Kum 39.43 

e(%) 80.81 

e(%) (%) for a time 5 times greater 35.67 
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Figure 5. Radius in pc of the PN Ring nebula as a function of the position angle in degrees. 
Physical parameters as in Table [T] 




(degree) 



Figure 6. Velocity in ^ of the PN Ring nebula as a function of the position angle in 
degrees. Physical parameters as in Table[T] 
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Figure 7. Map of the expansion velocity in ^ relative to the simulation of the PN Ring 
nebula when 300000 random points are selected on the surface. Physical parameters as in 
Table [J 



Table 4. Data of the simulation of the PN MyCn 18 



Initial expansion velocity ,Ro,kms [km s ^] 200 
Age(t4-to,4)[10^yr] 0.2 
Initial radius i?o,pc [pc] 0.001 
scaling h [pc] 1.0 x Rq 
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Figure 8. Continuous three-dimensional surface of the PN MyCn 18 : the three Eule- 
rian angles characterizing the point of view are $=130 °, B=40 ° and ^'=5 °. Physical 
parameters as in Table |4] 
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Figure 10. Radius in pc of the PN MyCn 18 as a function of latitude from 0° to 60° ( dotted 
line) when the physical parameters are those of Table |4] The points with error bar (1/10 of 
the value) represent the data of Table 1 in Dayal et al. 2000. 
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Figure 11. Map of the expansion velocity in ^ relative to the simulation of the PN 
MyCn 18 when 300000 random points are selected on the surface. Physical parameters as 
in Table m 



The combined effect of spatial asymmetry and field of velocity are reported in Figure 

m 

The efficiency of our code in reproducing the spatial shape over 12 directions of MyCn 
18 as given by formula (35 ) is reported in Table [5] This Table also reports the efficiency in 
simulating the shape of the velocity. 
Figure 



12 reports our results as well those of Table 1 in 1 46 1. 



4.3. Simulation of r^-Carinae in an exponentially varying medium 

A typical set of parameters which allows the Homunculus nebula around r/-Carinae to be 
simulated in the presence of a medium whose density decreases exponentially is reported 



Table 5. Reliability of the spatial and velocity shape of the PN MyCn 18. 



radius velocity 
eobs{%) 90.66 57.68 
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Figure 12. Velocity in ^ of the PN MyCn 18 as a function of the latitude in degrees when 
the physical parameters are those of Table[4] dotted line. The points with error bar (1/10 of 
the value) represent the data of Table 1 in Dayal et al. 2000. 



in Table |6] Table [7] presents numbers concerning the quality of fit. 

The bipolar character of the Homunculus is shown in Figure [13] In order to better 



visualize the two lobes, Figure 14 and Figure 15 show the radius and velocity as a function 



of the angular position 6. The orientation of the observer is characterized by the three Euler 
angles (<I>, 0, see (GT); different Euler angles produce different observed shapes. 



The velocity field is shown in Figure 16 



The accuracy with which our code reproduces the spatial shape and the velocity field 
over 18 directions of the Homunculus nebula as given by formula ( [35] ) is reported in Table[7] 
From a careful analysis of Table[7]we can conclude that the spatial shape over 18 directions 
is well modeled by an exponential medium , Eobs = 85%. The overall efficiency of the 



field is smaller Cobs = 75%. We can therefore conclude that formula (35 1 which gives the 



efficiency over all the range of polar angles represents a better way to describe the results 



in respect to the efficiency in a single direction as given by formula (34i. 



4.4. Simulation of r^-Carinae for a power law medium 

For assumed parameters see Table [6] Table J8] reports the accuracy of radius and velocity in 
two directions. 



4.5. Simulation of a spherical SNR , SN 1993J 

According to the power solution as given by ([4]) and to the data used in Section [3rL] e=98.54 
%. 
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Table 6. Parameter values used to simulate the observations of the hybrid Homunculus/ry- 
Carinae nebula for a medium varying exponentially (first 4 values) or a power law (2nd set 
of 4 values) 



Initial expansion velocity, -Ro.i [km s ^] 8000 

Age iU - toA) [10^ yr] 0.0158 

Scaling h [pc] 0.0018 

Initial radius Rq [pc] 0.001 



Initial expansion velocity, i?o,i [km s ^] 40, 000 
Age (t4 - to,4) [10^ yr] 0.0158 
Initial radius Rq [pc] 0.0002 
Power law coefficient a 2.4 



Table 7. Agreement between observations and simulations of the hybrid Homunculus/r;- 
Carinae nebula, for an exponentially varying medium. 



radius velocity 

e(%) — polar direction 97 99 

e(%) — equatorial direction 87 19 
eobs{%) 85 75 



Table 8. Agreement between model for a power law medium and observations for the 
hybrid Homunculus/ry-Carinae nebula. 



radius velocity 
e(%) — polar direction 78 76 
e(%) — equatorial direction 6 2 
eobs{%) 79 73 
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Figure 13. Simulations lead to this picture of the hybrid Homunculus/r^-Carinae nebula 
for an exponentially varying medium. The orientation of the figure is characterized by the 
Euler angles , which are 4>=130°, 0=40° and ^'=-140°. Physical parameters as in Table|6] 
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Figure 14. Radius of the hybrid HomuncuIus/r/-Carinae nebula as a function of latitude 
for an exponentially varying medium (dotted line) and astronomical data with error bar. 
Physical parameters as in Table [6] 




40 60 

(degree) 



Figure 15. Velocity of the hybrid Homunculus/r/-Carinae nebula as a function of latitude 
for an exponentially varying medium (dotted line) and astronomical data with error bar 
Physical parameters as in Table [6] 
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Figure 16. Map of the expansion velocity for an exponentially varying medium relative to 
the hybrid Homunculus/?7-Carinae nebula. Physical parameters as in Table |6] 
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Figure 17. The shape of the expanding envelope modeled by 2500 faces particularized for 
SN 1006 . The three Eulerian angles characterizing the point of view are <I'=75 °, 0=90 ° 
and ^=75 °. 



4.6. Simulation of a asymmetric SNR , SN 1006 

According to the numerical code developed in [68] in the case of a Gaussian profile we have 
e = 94.9% in the polar direction and e = 92.5% in the equatorial direction. From a practical 
point of view, the range of the polar angle 9 (180° ) is divided into ng steps and the range of 
the azimuthal angle (p (360° ) into n,/, steps. This yields (ng +1) (n^ +1) directions of motion 
which can also be identified with the number of vertexes of the polyhedron representing the 
volume occupied by the explosion ; this polyhedron varies from a sphere to an irregular 
shape on the basis of the swept-up material in each direction . In the plots showing the 
expansion surface of the explosion, the number of vertexes {n0+l)-{n^ +1) and the number 



of the faces ng • n^; are specified, for example in Figure 17 ng=50 and n^=50. 



5. Diffusion 

The mathematical diffusion allows us to follow the number density of particles from high 
values (injection) to low values (absorption). We recall that the number density is expressed 
in unit volume ^^'^ Symbol C is used in the mathematical diffusion and the symbol n in an 
astrophysical context. The density p is obtained by multiplying n by the mass of hydrogen 
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, rriH , and by a multiplicative factor , /, which varies from 1.27 in [69] to 1.4 in ||33| 

P = fmnn . (36) 

The physical process that allows the particles to diffuse is hidden in the mathematical dif- 
fusion. In our case the physical process can be the random walk with a time step equal to 
the Larmor gyroradius. In the Monte Carlo diffusion the step-length of the random walk is 
generally taken as a fraction of the side of the considered box. Both mathematical diffu- 
sion and Monte Carlo diffusion use the concept of absorbing-boundary which is the spatial 
coordinate where the diffusion path terminates. 

In the following, 3D mathematical diffusion from a sphere and ID mathematical as well 
Monte Carlo diffusion in presence of drift and are considered. 



5.1. 3D diffusion from a spherical source 

Once the number density , C, and the diffusion coefficient ,D, are introduced , Fick' s first 
equation changes expression on the basis of the adopted environment , see for example 
equation (2.5) in [ [70| . In three dimensions it is 

DV'^C , (37) 



dC 

'at 



where t is the time and is the Laplacian differential operator. 
In presence of the steady state condition: 

DV^C = . (38) 

The number density rises from at r=a to a maximum value Cm at r=b and then falls 
again to at r=c . The solution to equation p8|) is 



C(r) = A + 



B 



where A and B are determined by the boundary conditions , 



and 



Cab{r) = Cm[l 



a < r < b 



b < r < c 



(39) 



(40) 



(41) 



These solutions can be found in [70] or in 1 71 1 . 



5.2. ID diffusion with drift, mathematical diffusion 

In one dimension and in the presence of a drift velocity ,u, along the radial direction the 
diffusion is governed by Pick's second equation , see equation (4.5) in |70| , 



dC 



D 



Q^2 



dr 



(42) 
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Figure 18. Number density of the PN A39 as a function of the distance in pc from the 
injection when u = 1 , Cm = I, a = 69.6 arcsec, b = 87 arcsec , c = 89.6 arcsec and 
D = 2 (full line ), D = 7 (dashed ), D = 12 (dot-dash-dot-dash) and D = 17 (dotted ). 
The conversion from arcsec to pc is done assuming a distance of 2100 pc for A39. 

where u can take two directions. The number density rises from at r=a to a maximum 



value Cm at r=b and then falls again to at r=c . The general solution to equation (42) in 
presence of a steady state is 

C(r) = A + Bei" . (43) 

We now assume that u and r do not have the same direction and therefore u is negative ; 
the solution is 

C{r) = A + Be-T''' , (44) 

and now the velocity u is a scalar. 
The boundary-conditions give 

e — e D 

Ca,b,drift(.r) = C"™^^- txt: a <r <b downstream side , (45) 

e o — e o" 

and 

e D — e o 

Cb,c,drift{r) = Cm^^^ 6 < r < c upstream side . (46) 

e o — e 

A typical plot of the number density for different values of the diffusion coefficient is re- 



ported in Figure 1 8 



5.3. ID diffusion with drift, random walk 

Given a ID segment of length side we can implement the random walk with step-length A 

ndt 
X 



by introducing the numerical parameter NDIM = . We now report the adopted rules 



when the injection is in the middle of the grid : 
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1 . The first of the NPART particles is chosen. 

2. The random walk of a particle starts in the middle of the grid. The probabilities of 
having one step are pi in the negative direction (downstream) ,pi = ^ — /U x 5, and 
P2 in the positive direction (upstream) , P2 = \ + ^^ x \, where ^ is a parameter that 
characterizes the asymmetry (0 < ^ < 1). 

3. When the particle reaches one of the two absorbing points , the motion starts another 
time from (ii) with a different diffusing pattern. 

4. The number of visits is recorded on 7W , a one-dimensional grid. 

5. The random walk terminates when all the NPART particles are processed. 

6. For the sake of normalization the one-dimensional visitation or number density grid 
M is divided by NPART. 

There is a systematic change of the average particle position along the x-direction: 

{dx) = /i A , (47) 

A 

Vtr 

asymmetry ,/i , that characterizes the random walk is 



for each time step. If the time step is = where vtr is the transport velocity, the 



/i = — . (48) 

Vtr 



Figure 19 reports M.{x), the number of visits generated by the Monte Carlo simulation as 



well as the mathematical solution represented by formulas ( |45| ) and (46 1 



The solutions of the mathematical diffusion equations (45 1 and (46 1 can be rewritten at 
the light of the random walk and are 

e A " — e A ' 

Ca,b,Mc{f) = Cm — 271 27]— a <r <b downstream side , (49) 



and 



e \°- — e >^ 



e A — e A ' 



c,Mc{f) — Cm _2£, _2m, b <r < c upstream side . (50) 

e~~^ — e~~ 

6. Radiative transfer equation 

The transfer equation in the presence of emission only , see for example ||72j or | [73| , is 

^ = -kXL+j.C , (51) 

ds 

where /j^ is the specific intensity , s is the line of sight , j,^ the emission coefficient, a 
mass absorption coefficient, ( the mass density at position s and the index v denotes the 



interested frequency of emission. The solution to equation ( 5 1 1 is 



/,(r,) = ^(l-e---^(^)) , (52) 




Figure 19. Number density relative to PN A39 of the ID asymmetric random walk (full 
line), NDIM=401 ,NPART=200 ,side = 40 arcsec , A = 0.1 arcsec and ^ =- 0.013. For 
astrophysical purposes n is negative. The theoretical number density as represented by 
formulas (45 1 and (46 1 is reported when u = 1 , Cm = 1, a = 60 arcsec, 6 = 80 arcsec 
, c = 100 arcsec and D = 3.84 (dotted line ). The conversion from arcsec to pc is done 
assuming a distance of 2100 pc for A39. 
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where Ty is the optical depth at frequency u 

dTy = kyC,ds . (53) 

We now continue analyzing the case of an optically thin layer in which t,^ is very small ( or 
ky very small ) and the density is substituted with our number density C(s) of particles. 
Two cases are taken into account : the emissivity is proportional to the number density and 
the emissivity is proportional to the square of the number density . In the linear case 

j.( = KC{s) , (54) 

where K is a constant function. This can be the case of synchrotron radiation in presence 
of a isotropic distribution of electrons with a power law distribution in energy, N{E), 

N{E)dE = K,E-'<f , (55) 

where Kg is a constant. In this case the emissivity is 

j> pa 0.933 X 10-23a(7^)i^,i7 2 ) 2 ^ (55) 

where v is the frequency and 0(7^) is a slowly varying function of 7^ which is of the order 
of unity and is given by 

a(7/) = 2fa-3)/^ ^/ + V\ ipi^)ripl±l) , (57) 

yu) 7J + 1 ^ 12 ^ ^ 12 ^ ' ^ ' 



for 7j > 2> see formula (1.175 ) in |74| . The synchrotron emission is widely used to 
explain the radiation observed in SNR, see |75 - 80 1. This non thermal radiation continuum 
emission was also detected in a PN associated with a very long-period OH/IR variable star 
(V1018 Sco), see [81 1. 
In the quadratic case 

juC = K2C{sf , (58) 
where K2 is a constant function. This is true for 



• Free-free radiation from a thermal plasma, see formula ( 1 .2 1 9) in 1 74 1 . This radiation 
process was adopted by |[48| in the little Homunculus. 



• Thermal bremsstrahlung and recombination radiation , see formula (1.237) in [ 74] 
This radiation process was adopted in PNs by [[82H84). 



The intensity is now 

Iu{s) = K / C{sf)dsf optically thin layer linear case , (59) 

J Sn 



or 



Iu{s) = K2 I C{sf)'^dsf optically thin layer quadratic case 

J so 



(60) 
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In the Monte Carlo experiments the number density is memorized on the grid A4 and the 
intensity is 

-^(^)i) = ^ ^ ■s X -^(^) J) ^) optically thin layer linear case , (61) 

k 

or 

Hhj) = ''^^sxM.{i,j,k)'^ optically thin layer quadratic case , (62) 

k 

where As is the spatial interval between the various values and the sum is performed over 
the interval of existence of the index k. The theoretical intensity is then obtained by inte- 
grating the intensity at a given frequency over the solid angle of the source. 

In order to deal with the transition to the optically thick case, the intensity is given by 

Hhj) = -^(1 -exp(-i^,^ As X S{iJ,k))) (63) 

" k 

Thin I — Thick , 

where Ka is a constant that represents the absorption. Considering the Taylor expansion of 
the last formula (63 1, equation (61 1 is obtained. 



7. Images 

The image of a PN , r/-Carinae and a SNR can be modeled once an analytical or numerical 
law for the intensity of emission as a function of the radial distance from the center is given. 
Simple analytical results for the radial intensity can be deduced in the rim model when the 
length of the layer and the number density are constants. The integration of the solutions to 
the mathematical diffusion along the line of sight allows us to deduce analytical formulas 
in the spherical case. The complexity of the intensity in the aspherical case can be attached 
only from a numerical point of view. 



7.1. 3D Constant Number density in a rim model 

We assume that the number density C is constant and in particular rises from at r = a 
to a maximum value Cm , remains constant up to r = 6 and then falls again to 0. This 



geometrical description is reported in Figure 20 The length of sight , when the observer 



is situated at the infinity of the x-axis , is the locus parallel to the x-axis which crosses the 
position y in a Cartesian x — y plane and terminates at the external circle of radius b. The 
locus length is 

loa = 2x (y^b^ - y2 _ _ ^2) ;0<y<a 

/,fe = 2x (V&2_y2) .a<y<b . (64) 

When the number density Cm is constant between two spheres of radius a and b the intensity 
of radiation is 

/Oa = Cm X 2 X ( _ y2 _ ^^2 _ ^2) ;0<y<a 

lab = Cm X 2 X iVb^-y^) ■,a<y<b . (65) 
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Figure 20. The two circles (section of spheres) which include the region with constant 
density are represented through a full Une. The observer is situated along the x direction, 
and three Unes of sight are indicated. 
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1 ■ 1 1 1 1 1 r 




Arcsec from Nebula Center 



Figure 21. Cut of the mathematical intensity / of the rim model ( equation (65 1) crossing 
the center (full line ) of the PN A39 and real data (dotted line with some error bar ) . The 
number of data is 801 and for this model = 1.487 against = 0.862 of the rim model 
fully described in Jacoby et al. (2001). 



The comparison of observed data of A39 and the theoretical intensity is reported in Fig- 
ure |2T] when data from Table |9] are used. 

The ratio between the theoretical intensity at the maximum , {y = b) , and at the 
minimum , (y = 0) , is given by 

i{y = b) .... 



7.2. 3D diffusion from a sphere, square dependence 



Figure 22 shows a spherical shell source of radius b between a spherical absorber of radius 
a and a spherical absorber of radius c. 

The number density rises from at r=a to a maximum value Cm at r=b and then falls 
again to at r=c . 



The numbers density to be used are formulas (40 1 and (41 



imposed ; these two numbers density are inserted in formula ( [58] ) which represents the 
transfer equation with a quadratic dependence on the number density. An analogous case 



once r 



y2 is 



was solved in 1 85 1 by adopting a linear dependence on the number density . The geometry 
of the phenomena fixes three different zones (0 — a, a — b, b — c) in the variable y, ; the first 
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Table 9. Simulation of the PN A39 with the rim model 



symbol 


meaning 


value 


a 


radius of the internal sphere 


72.5" 


b 


radius of the external sphere 


90.18" 


Rshell 


observed radius of the shell 


77" 


S r shell, t 


theoretical thickness of the shell 


17.6" 


S rshell 


observed thickness of the shell 


10.1" 


^center 


ratio of observed intensities 


(1.88 - 2.62) 


^max 

l(y=0) 


ratio of theoretical intensities 


3.03 



I 



iC=0 




Figure 22. The spherical source inserted in the great box is represented through a dashed 
line, and the two absorbing boundaries with a full line. The observer is situated along the x 
direction, and three lines of sight are indicated. Adapted from Figure 3.1 by Berg (1993) . 



piece , l\y) , is 
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I\y) = / , 2C„\dx+ / 2Cldx 

-2ayln(V'a2-i/2 + a)6^ + 2a^ arctanf-^^ —)ch 

y 

+2 yy/W^cb + 2 a ln( V?'^ - y2 + b)yb'^ + 2 a ln( V?'^ _ y2 + ^^^^^ 
-2cyln(V?'2-y2 + _ 2c2 arctan(^^?^?)6a - 2y^h'^-y'^ba - 

y 

2cyln{^/b'^ - y'^ + h)a^ + 2cln(v^c2 - y2 + c)yb^ 

+2cln(vc2^-y2 + c)y(]? + 2 \J — y^yba + 2(? arctan(-^^ )ba 

—2 ay ln(-\/ a? — y2 -|- 0)^2 _|_ ^ ^2 _ y'^y(^ 
— 4cln(i/ c2 — j/2 -|- c)y6a + Aay ln(\/ a2 — y2 _|_ (j)gj, 

- V c2 — y^ya^ + arctan( — y v 62 — ?/2g2 _ ^^2 

_ y2 

— \J(? — y'^yh^ + y\Jb'^ — y'^a^ + arctanf^^^ 

y 

y 



.c2 arctan(^^^^)62 + arctan( ^^ ^^ )c^1 (67) 

y y 

Q<y <a 
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The second piece , (y) , is 



f\y) = / 2C„Va;+ / 2Cldx 



2 _ „,2 



-c^ arctan( ^ 



-yV^^^a' + arctan(^^?^)a2 + arctan( ^''' + y^?^^.,.^ 



c — y ya 

y ' ' y ' 



2 



+ a/c2 - + 2 a ln(y)2/6^ + 2 a ln(y)yc^ + 2 cy ln{^/b'^ - y"^ + b)l) 

-2 arctan(^^^?^5^)c6 - 2 y^ly^ - y^cb -2a ln( V^P^^ + b)yb^ 



— 2aln(vfc^ — + ^'jyc^ — 2c^arctan( )ba — 2\J(? — y^yba 



/■U2 _ ..2 

-2cln(\/c2 -y2 + c)yo^ - 2cln(Vc2 -y^ + c)yb^ + 2 arctan( ^ ^ ^ )ba 

+2y^/W^^ba + 2cyhi{^/b'^ - y'^ + 6)a^ -Aahi{y)ycb + Ac\u{^/c^ - + c)y6a] (68) 

a < y < b . 
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The third piece , P^\y) , is 



Jo 



y(62-25a + a2)(c2-2c6 + 62) 

— V o2 — y'^yb^ + (? arctan( 

V 



+0? arctan( ^ ^ )fe2 + yjb"^ — y'^yc? — yjb"^ — y'^ya^ — (? arctan( 

y y 





y2 


y 








y2 



-a^ arctan( — arctan( + 2 cy \n{\/b'^ — + b)a^ 



y y 



+2 ay ln(\/a2 - y^ + a)c^ + 2 arctan(-^^ —)cb - 2a\n{^Jb'^ - y^ + b)yc^ 



y 



-2 a ln(V62 _ y2 + ^)y52 _ 2 ^52 _ ^2^^^ _ 2 arctan(^^^^^' — ^)c6 
+2cyln(\/62 _ y2 _^ 5)52 _ 2cln(v^c2 - y2 + c)y6^ - 2cln(v^c2 - y'^ + c)ya^ 



,2 _ „2 



-2 c2 — y'^ba — 2(? arctanf^^^- —)ba + 2 ay ln(\/ a? — y2 + a}^"^ 

y 



.2 _ ,,2 



fv_ _ , _ 

+2 \/ — y^ycb + 2 \/b'^ — y'^yba + 2 arctan( )6a 



,2 _ y2 



+4cln(A/c2 - 2/2 + c)y6a - 4 ay 111(1/02 - y2 + a)c6] (69) 

b < y < c . 



The profile of / made by the three pieces ([68]), ( [69| ) and ([70|), can be calibrated on the 
real data of A39 and an acceptable match is realized adopting the parameters reported in 
Table [lOl 

The theoretical intensity can therefore be plotted as a function of the distance from the 



center , see Figure 23 or as an image , see Figure 24 

The effect of the insertion of a threshold intensity , Itr, given by the observational 
techniques , is now analyzed. The threshold intensity can be parametrized to Imax^ the 
maximum value of intensity characterizing the ring: a typical image with a hole is visible 



in Figure |25j when Itr = Imax/2. 

The position of the minimum of / is at y = and the position of the maximum is 
situated at y = 6. 

The ratio between the theoretical intensity at maximum , Imax y = b , and at the 
minimum (y = 0) is given by 

Imnx Numerator ^^^^ 



I{y = 0) Denominator 
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Table 10. Simulation of the PN A39 with 3D diffusion 



symbol 


nieamng 


value 


a 


radius of the internal absorbing sphere 


65.96" 


b 


radius of the shock 


80" 


c 


radius of the external absorbing sphere 


103.5" 


Rshell 


observed radius of the shell 


77" 


S r shell 


observed thickness of the shell 


10.1" 


^center 


ratio of observed intensities 


(1.88 - 2.62) 


Imax 

Hy=o) 


ratio of theoretical intensities 


2.84 




Arcsec from Nebula Center 



Figure 23. Cut of the mathematical intensity / (formulas ( 68 1, ( 69 1 and ( 70l) , crossing 
the center (full line ) of PN A39 and real data (dotted line with some error bar). The 
number of data is 801 and for this model = 19.03 against = 12.60 of the rim model 
fully described in Jacoby et al. (2001). 
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Figure 25. The same as Figure 



24 



but with Itr = /max/2 
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where 



Numerator = (b^ — 2ba + 



; (2 cb ln(6) -2c ln{\/ - 5^ + c)b + 6\/ - 62 + c2 a,rctani 



62 



)) 



(71) 



and 



Denominator = (72) 
2 6(o^c -c^a-2 bca ln(a) + 2 bca ln(c) - 6a^ + 6c^ - 6^c + 6^a + 6^a ln(a) 
-c^a ln(6) + 62cln(6) - 6^aln(6) - 62cln(c) + a^c ln(6) + c2aln(a) - a'^cln{c)) . 



The ratio rim(maximum) /center(minimum) of the observed intensities as well as the 



theoretical one are reported in Table 10 for A39 |l86j 



7.3. 3D diffusion from a sphere, linear dependence 

The concentration rises from at r=a to a maximum value Cm at r=b and then falls again 
to at r=c. The concentrations to be used are formulas ( 40 1 and ( 4 1 1 once r = a/x2 + y2 



imposed; these two concentrations are inserted in formula ( 52 1 which represents the transfer 
equation. The geometry of the phenomenon fixes three different zones (0 — a, a — 6, 6 — c) 
for the variable y, see |85 86J; the first segment, (y), is 



^ bCm\/o? - 
-b + a 



bCma In ( \/a2 - y2 _^ 



bCmVb^ 



-b + a 



-b + a 



bCmaln(./W^ + b) bC^c\u(^/W^ + b 
+2 ^ ^ + 2 



-h + a 



-c + b 



_^bCrnV^ZZ _ ^ ^CraCln[y^C^-y^ + c) ^ ^^y^^/^—y2 



-C + b 



-c + b 



-c + b 
<y < a 



(73) 



The second segment, I (y), is 



bCmaln (y^ 
-b + a 



bCmVb^ 



y 



-b + a 



bCma\n(.JW^ + b] bC^c\u(.JW^ + b 

+2 ^ ^ + 2 



-b + a 



-c + b 



hCraVb^ 



-c + b 



6C„cln(Vc2-y2 + cj f^—^ 
2 ^ ^ + 2 — ^ , , ^ (74) 



-c + b 



-c + b 
a < y < b 
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Table 11. Simulation of the SNR SN 1993 J by 3D diffusion, optically thin case 

symbol meaning value 

a radius internal sphere 1.76{mas) 

b radius shock 2.2{mas) 

c radius external sphere 5.0{mas) 

ratio observed intensities 1.7926 



^center 
I max 

l{y=0) 



ratio theoretical intensities 1.7927 



The third segment, I^^^{y), is 



y = - 2 ^— ^ + 2 

— c + —c + b 



bCmV^ 

-c + b 
b < y < c 



y 



(75) 



The profile of / made up of the three segments (73 1, (74i and (76 1, can be calibrated 
against the real data of SN 1993J and an acceptable match can be achieved by adopting 



the parameters reported in Table 1 1 The theoretical intensity can therefore be plotted as a 



function of the distance from the center, see Figure 26 or as a contour map, see Figure 27 



The position of the minimum of / is at y = and the position of the maximum is 
situated in the region a < y < b, or more precisely at: 



{b — 2a + c) a {ab — 2bc + ac) 



b-2a+c 



(76) 



This means that the maximum emission is not at the position of the shock, identified here 



as 6, but shifted a little towards the center; see Figure 28 



The ratio between the theoretical intensity at maximum , I max , as given by formula ( 76 1 
and at minimum (y = 0) is given by 



In 



Numerator 



I{y = 0) Denominator 



(77) 



42 



Zaninetti 




Radius (mas) 



Figure 26. Cross-section of the mathematical intensity / (formulas (73 1, (74i and (76 1), 
through the center (dotted line) of SN 1993J and real data (empty stars), = 100.49 The 
real data made on day 1889 of SNR SN 1993J after the explosion have been extracted by 



the author from Figure 3 of Marcaide et al. (2009). Parameters as in Table 1 1 
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Figure 27. Contour map of / adjusted to simulate the SNR SN 1993 J . Parameters as in 
Table [m 
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Figure 28. Cross-section through the mathematical intensity / towards the edge of the 
SNR SN 1993J . The three parameters which characterize the expanding PN, a , b and c, 
are reported. Parameters as in Table [TT] 



where 



Numerator = 

, , (ac + ab — 2 cb) a\ _ / (ac + ab — 2 cb) a\ , ^ (c + b)(b — a)'^ 
-ami — - — ; — ] c + am { — - — ; — b — 2 \ — — —c 



b+c-2a J V b+c-2a J V b+c-2a 



2cln I J(^±^H^ + , K + 2 . M^,. 
b+c-2a I V b+c-2a 



, , (a-c?(c + b) 



, , (a-cy(c + b) \ (a-cy(c + b). 



V b + c — 2a 
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Figure 29. Cross-section through the mathematical intensity / (formulas (73 1, (74 1 and 
(76 1), in the optically thin case (dashed line, = 237.3), and optically thick case (full 
line, = 84.7) and real data (empty stars) for SNR SN 1993J . Parameters as in Table 12 



and 



Denominator 

-2acln (a) + 2 6a In (a) - 2 6a In (6) + 2 6cln(6) - 2 6c In (c) + 2acln (c) 



(79) 



The observed ratio as well as the theoretical ratio are reported in Table 1 1 



The effect of absorption is easily evaluated by applying formula ( 63 1 and fixing the 



value of Ka- The result is shown in Figure 29 



7.4. 3D complex morphologies of PN 



The numerical approach to the intensity map can be implemented when the ellipsoid that 
characterizes the expansion surface of the PN has a constant thickness expressed , for ex- 
ample , as rmin/ f where rmin is the minimum radius of the ellipsoid and / an integer. 
We remember that / = 12 has a physical basis in the symmetrical case , see [33J . The 
numerical algorithm that allows us to build the image is now outlined 

• A memory grid A^(i, j, k) that contains NDIM^ pixels is considered 

• The points of the thick ellipsoid are memorized on M 

• Each point of Jv[ has spatial coordinates x,y,z which can be represented by the 
following 1x3 matrix ,A, 

X 



A 



y 



(80) 
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Table 12. Simulation of the SNR SN 1993J with 3D diffusion, optically thick case with 
Ka = 0.2. 



symbol 


meaning 


value 


a 


radius internal sphere 


2.0l{mas) 


b 


radius of shock 


2.2{mas) 


c 


radius external sphere 


5.0{mas) 


I center 


ratio observed intensities 


1.7926 


I max 

l(y=0) 


ratio optically thin case 


2.0491 


I max 

l{y=0) 


ratio optically thick case 


1.6741 



The point of view of the observer is characterized by the Eulerian angles ($, 0, 'J/) 
and therefore by a total rotation 3x3 matrix , E , see \61} . The matrix point is now 
represented by the following 1x3 matrix , B, 



B = E-A 



(81) 



• The map in intensity is obtained by summing the points of the rotated images along 
a direction , for example along z , ( sum over the range of one index, for example k ). 



Figure 30 reports the rotated image of the Ring nebula and Figure 3 1 reports two cuts along 



the polar and equatorial directions. 



Figure 32 reports the comparison between a theoretical and observed east-west cut in 
Hjs that cross the center of the nebula, see Figure 1 in [87| |. 

A comparison can be made with the color composite image of Doppler-shifted H2 emis- 
sion as represented in Figure 2 in ||42|. 

In order to explain some of the morphologies which characterize the PN's we first map 



MyCn 1 8 with the polar axis in the vertical direction , see map in intensity in Figure 33 The 



vertical and horizontal cut in intensity are reported in Figure 35 The point of view of the 
observer as modeled by the Euler angles increases the complexity of the shapes : Figure 34 



reports the after rotation image and Figure 36 the vertical and horizontal rotated cut. The 
after rotation image contains the double ring and an enhancement in intensity of the central 
region which characterize MyCn 18. 

This central enhancement can be considered one of the various morphologies that the 
PNs present and is similar to model BLi — F in Figure 3 of the Atlas of synthetic line 
profiles by p6|. 



7.5. 3D complex morphology of the hybrid r^-Carinae 



Here we adopt the numerical algorithm developed in the previous Section 7.4. An ideal 



image of the Homunculus nebula having the polar axis aligned with the z-direction which 



means polar axis along the z-direction, is shown in Figure 37 and this should be compared 
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-0.1 0.1 



arcsec from center 

Figure 30. Map of the theoretical intensity of the PN Ring nebula. Physical parameters as 
in Table [T] and /=12 . The three Eulerian angles characterizing the point of view are <I>=180 
°, e=90 ° and ^=-30 °. 




-OJ -0.t 0.1 0^ 

R (pc) 

Figure 31. Two cut of the mathematical intensity / crossing the center of the PN Ring 
nebula: equatorial cut (full line) and polar cut (dotted line) . 
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Figure 32. Cut of the mathematical intensity / of the PN Ring Nebula crossing the center 
(full line ) and real data of Hp (dotted line with some error bar ) . The number of data is 
250 and for this model = 15.53 . The real data are extracted by the author from Figure 1 
of Gamett and Dinerstein 2001. 



with the H2 emission structure reported in Figure 4 of |[55|. A model for a realistically 



rotated Homunculus is shown in Figure 38 This should be compared with Figure 1 in 
or Figure 1 in | ,53| . 

The rotated image exhibits a double ring and an intensity enhancement in the central 



region which characterizes the little Homunculus, see | |47|[48|pl|[89| . Figure 39 and Fig- 
ure [40] show two cuts through the Homunculus nebula without and with rotation. The 



intensity enhancement is due to a projection effect and is an alternative for the theory that 
associates the little Homunculus with an eruption occurring some time after the Great Erup- 
tion, see pT]|89[ . We briefly recall that a central enhancement is visible in one of the vari- 
ous morphologies characterizing planetary nebulae. This can be compared with the model 
BLi — F in Figure 3 of the Atlas of synthetic line profiles by |fT6l. 

Such cuts are common when analyzing planetary nebulae. As an example Figure 4 
in | |40| reports a nearly symmetrical profile of the intensity in the [OIII] image of A39, 
a nearly spherical planetary nebula. Another example is the east-west cut in Hfi for the 
elliptical Ring nebula, crossing the center of the nebula, see Figure 1 in [87] . Such intensity 
cuts are not yet available for Ty-Carinae and therefore can represent a new target for the 
observers. 



7.6. 3D complex morphology of a SNR , SN 1006 

The theory of an asymmetric SNR was developed in Sect. 4. 1 of [68 1 in which an expansion 
surface as a function of a non-homogeneous ISM was computed: in the same paper Figure 8 
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'-0.2 -0.1 0.1 0.2 



crcsec from center 

Figure 33. Map of the theoretical intensity of the PN MyCn 18 . Physical parameters as in 
Table [6] and /=12 . The three Eulerian angles characterizing the point of view are ^>=180 ° 
, e=90 ° and ^-=0 °. 



Interaction of PN and SNR with the ISM 



49 




-0.1 0.1 

arcsec from center 



Figure 34. Map of the theoretical intensity of the rotated PN MyCn 1 8 . Physical parame- 
ters as in Table |6] and /=12 . The three Eulerian angles characterizing the point of view are 
^>=130 °, e=40 ° and ^=5 °. 




R (pc) 



Figure 35. Two cut of the mathematical intensity / crossing the center of the PN MyCn 18 
: equatorial cut (full line) and polar cut (dotted line) . Parameters as in Figure [33 
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Figure 36. Two cut of the mathematical intensity / crossing the center of the rotated PN 
MyCn 18 nebula: equatorial cut (full line) and polar cut (dotted line) . Parameters as in 



Figure 34 



models SN1006. The diffusing algorithm adopted here is the 3D random walk from many 
injection points ( in the following IP) 

1 . The first IP is chosen 

2. The first of the NPART electrons is chosen. 

3. The random walk of an electron starts where the selected IP is situated. The electron 
moves in one of the six possible directions. 

4. After N steps the process restarts from (2) 

5. The number of visits is recorded on Ai"^ , a three-dimensional grid. 

6. The random walk terminates when all the NPART electrons are processed. 

7. The process restarts from (1) selecting another IP 

8. For the sake of normalization the one-dimensional visitation/concentration grid 
is divided by NPART. 

The IP are randomly selected in space, and the radius is computed by using the method 
of bilinear interpolation on the four grid points that surround the selected latitude and lon- 
gitude, ( 1 66 1 ). The radius will be the selected value + R/24 in order to generate the IP 
where the action of the shock is maximum. 

Our model gives radial velocities , Vtheo > 2211 km s^^ < Vtheo ^ 3580 km 



and the map of the expansion velocity is reported in Figure 41 from which it is possible to 
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arcsec from cente 



Figure 37. Map of the theoretical intensity of the hybrid Homunculus/ry-Carinae nebula 
in the presence of an exponentially varying medium. Physical parameters as in Table |6] 
The three Euler angles characterizing the orientation are <I>=180°, 0=90° and ^'=0°. This 
combination of Euler angles corresponds to the rotated image with the polar axis along the 
z-axis. 
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Figure 38. Model map of the of the hybrid Homunculus/?7-Carinae nebula rotated in accor- 
dance with the observations, for an exponentially varying medium. Physical parameters as 
in Table[6] The three Euler angles characterizing the orientation of the observer are <^>=130°, 
0=40° and ^'=-140°. This combination of Euler angles corresponds to the observed image. 
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Figure 39. Two cuts of the model intensity across the center of the hybrid Homunculus/77- 
Carinae nebula for an exponentially varying medium: equatorial cut (full line) and polar cut 
(dotted line). Parameters as in Figure 37 
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Figure 40. Two cuts of the model intensity across the center of the realistically rotated 
hybrid Homunculus/r/ — car nebula for an exponentially varying medium: equatorial cut 



(full line) and polar cut (dotted line). Parameters as in Figure 38 



visualize the differences in the expansion velocities among the various regions as well as 
the overall elliptical shape. 

Before continuing we should recall that in the presence of discrete time steps on a 3D 



lattice the average square radius ,{R'^{N)), after N steps (see |90|, equation (12.5 )) is 

{R^{N)) ~ 6DN , (82) 
from which the diffusion coefficient , D ,is derived 

The two boundaries in which the random walk is taking place are now represented by two 
irregular surfaces. It is possible to simulate them by stopping the random walk after a 
number of iterations N given by 

N = NINT(^\f , (84) 

^ 24: 6 

where Rpc represents the averaged radius in pc . These are the iterations after which ac- 



cording, to formula (83 1, the walkers reach the boundaries at a radial distance given by 

from the place of injection; in other words we are working on an unbounded lattice 
. The influence of velocity on the flux F of radiation can be inferred from the suspected 
dependence when non-thermal emission is considered, see equation (9.29) in | [9T| , 

1 o 

F = Xt^f^HnoV^ , (85) 
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Figure 41. Map of the expansion velocity relative to the simulation of SNR SN 1006 when 
190000 random points are selected on the surface. The physical parameters aie the same as 
in Figure 8 of 1 ,68J . 



where xt represents the efficiency of conversion of the unitarian flux of kinetic energy, 
the mass of the hydrogen nucleus, no the particles/cm^ and Vs the velocity of the shock. 

Assuming that the flux reversed in the non-thermal emission follows a similar law 
through the parameter xx ( the efficiency in the X-region) the effect of velocity is sim- 
ulated through the following algorithm. Once the IP are spatially generated, the number of 
times NTIMES over which to repeat the cycle is given by 

NTIMES = 1 + NTIMESmax * ( " ~ ^"'^ f , (86) 

where NTIMESmax is the maximum of the allowed values of NTIMES minus 1 , and v 
is the velocity associated to each IP. The asymmetric contour map obtained when the spatial 



step is PS 2 * gyro-radius is reported in Figure 42 and the cut along two perpendicular lines 



of the projection grid in Figure 43 



In Figure 43 the asymmetry both in the peak to peak distance and the difference in the 
two maximum is evident. The ratio between the X-ray emission in the bright limbs (NE or 
SW) and toward the northwest or southeast ( at 2 keV) is around 10 , see Figure 5 top right 



in | [92| . Conversely our theoretical ratio , see Figure|43j is 9.84. It is also possible to plot the 
maximum of the theoretical intensity as function of the position angle , see Figure 44 The 
reader can make a comparison with the observational counterpart represented by Figure 5 
top right, dashed line in [92] . 



8. Conclusions 

Law of motion The law of motion in the case of a symmetric motion can be modeled by a 
power law solution , the Sedov Solution or the radial momentum conservation. These three 
models allow to determine the approximate age of A39 which is 8710 yr for the Sedov 
solution and 50000 yr for the radial momentum conservation. In presence of gradients 
as given , for example , by an exponential behavior, the solution is deduced through the 
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pc 



Figure 42. Contour of the intensity / in the X-rays of SNR SN 1006 . The parameters are 
sidesArR=18.37 pc , S = 6.12 IQ-^ pc, p = 2.8 IQ-^ pc, NDIM=3001, IP = 190000^ , 
NPART= 100, NTIMESmax='^^, great box. Optically thin layer. 




pc 



Figure 43. Two cut along perpendicular lines of / for SNR SN 1006 reported in Figure 
Optically thin layer. 
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Figure 44. Azimuthal maximum of the profiles of intensity as function of the position 
angle in degrees for SNR SN 1006 . Same parameters as in Figure 42 
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radial momentum conservation. The comparison with the astronomical data is now more 
complicated and the single and multiple efficiency in the radius determination have been 
introduced. When , for example , MyCn 18 is considered , the multiple efficiency over 
18 directions is 90.66% when the age of 2000 yr is adopted. In the case of ry-Carinae the 
multiple efficiency over 18 directions is 85% for a fixed age of 158 yr in an exponentially 
varying medium. In the case of the weakly asymmetric SN 1006 the efficiency is 94.9% 
in the polar direction and 92.5% in the equatorial direction for a fixed age of 974 yr in an 
exponentially varying medium. 
Diffusion 

The number density in a thick layer surrounding the ellipsoid of expansion can be con- 
sidered constant or variable from a maximum value to a minimum value with the growing 
or diminishing radius in respect to the expansion position. In the case of a variable number 



density the framework of the mathematical diffusion has been adopted, see formulas (40 1 



and (41 1. The case of diffusion with drift has been analytically solved , see formulas (45 1 
and ( 46 1, and the theoretical formulas have been compared with values generated by Monte 
Carlo simulations. 

Formation of tlie image 

The intensity of the image of a symmetrical PN or SNR in the case of optically thin 
medium can be computed through an analytical evaluation of lines of sight when the num- 



ber density is constant between two spheres see formula (65 1. The case of a symmetrical 



diffusive process which is built in presence of three spheres we should distinguish between 
• intensity of emission proportional to the square of the number density corresponding 



to the case of PN , see formulas ( 69 1 and ( 70 1. 



intensity of emission proportional to the number density corresponding to the case of 
SNR, see formulas (173]), {U^ and (l76ll. 



In the case of complex morphologies assuming an optically thin medium, it is possible 
to make a model image once two hypotheses are made: 

1. The thickness of the emitting layer, Ai?, is the same everywhere Ai? = 0.03i?max> 
where Rmax is the maximum radius of expansion. 

2. The density of the emitting layer is constant everywhere 

A 2D image of the PNs Ring nebula and MyCn 1 8 , the hybrid Homunculus/ry-Carinae neb- 
ula and the weakly asymmetric SNR SN 1006 are shown respectively in Figures 30j 33 37 
and 021 
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